MCSCF
Innehåll
MCSCF¶
The core functionality of MultiPsi is multiconfigurational methods, with multiconfigurational self-consistent field (MCSCF) as its cornerstone.
CASSCF¶
The most commonly used form of MCSCF nowadays is the complete active space self-consistent field (CASSCF), which is a form of MCSCF where the configurations are generated from a full CI expansion within an active space consisting of a few orbitals. Because of the factorial growth of full CI with the number of active orbitals, the active space should usually be limited to less than 20 electrons in 20 orbitals, with the current record being a CAS(22,22) on the triplet state of hydrogenase. However, this large calculation was only possible on a large cluster using the entire distributed memory available.
To define a CASSCF, we simply need to define our OrbSpace as a full CI expansion (see previous chapters) and provide it to our MCSCF module. For example, we can start with a simple Hartree-Fock, here for furan:
import veloxchem as vlx
import multipsi as mtp
furan_xyz="""9
C -0.86213 -0.90784 0.00007
H -1.63433 -1.64264 -0.00003
C 0.50727 -0.90524 0.00007
C 0.92057 0.47886 -0.00003
C -0.22323 1.23186 -0.00003
O -1.35123 0.40376 -0.00013
H 1.17117 -1.74724 0.00017
H 1.93767 0.81866 0.00007
H -0.46573 2.26986 -0.00013
"""
molecule = vlx.Molecule.from_xyz_string(furan_xyz)
basis = vlx.MolecularBasis.read(molecule,"def2-sv(p)")
scfdrv = vlx.ScfRestrictedDriver()
scfdrv.compute(molecule, basis)
* Warning * Environment variable OMP_NUM_THREADS not set.
* Warning * Setting OMP_NUM_THREADS to 16.
Self Consistent Field Driver Setup
====================================
Wave Function Model : Spin-Restricted Hartree-Fock
Initial Guess Model : Superposition of Atomic Densities
Convergence Accelerator : Two Level Direct Inversion of Iterative Subspace
Max. Number of Iterations : 50
Max. Number of Error Vectors : 10
Convergence Threshold : 1.0e-06
ERI Screening Scheme : Cauchy Schwarz + Density
ERI Screening Mode : Dynamic
ERI Screening Threshold : 1.0e-12
Linear Dependence Threshold : 1.0e-06
* Info * Nuclear repulsion energy: 158.8888494994 a.u.
* Info * Overlap matrix computed in 0.00 sec.
* Info * Kinetic energy matrix computed in 0.00 sec.
* Info * Nuclear potential matrix computed in 0.00 sec.
* Info * Orthogonalization matrix computed in 0.00 sec.
* Info * SAD initial guess computed in 0.00 sec.
* Info * Starting Reduced Basis SCF calculation...
* Info * ...done. SCF energy in reduced basis set: -228.337309732541 a.u. Time: 0.27 sec.
* Info * Overlap matrix computed in 0.00 sec.
* Info * Kinetic energy matrix computed in 0.00 sec.
* Info * Nuclear potential matrix computed in 0.00 sec.
* Info * Orthogonalization matrix computed in 0.00 sec.
Iter. | Hartree-Fock Energy | Energy Change | Gradient Norm | Max. Gradient | Density Change
--------------------------------------------------------------------------------------------
1 -228.430357327904 0.0000000000 0.16297779 0.01240917 0.00000000
2 -228.433229998567 -0.0028726707 0.05564260 0.00546548 0.06351427
3 -228.433553776745 -0.0003237782 0.00841002 0.00075289 0.01658680
4 -228.433566986231 -0.0000132095 0.00299022 0.00026049 0.00498463
5 -228.433568460937 -0.0000014747 0.00114921 0.00013437 0.00167819
6 -228.433568710661 -0.0000002497 0.00030252 0.00003062 0.00054130
7 -228.433568755136 -0.0000000445 0.00007367 0.00000963 0.00033342
8 -228.433568757937 -0.0000000028 0.00002138 0.00000334 0.00008594
9 -228.433568758029 -0.0000000001 0.00001803 0.00000167 0.00001459
10 -228.433568758067 -0.0000000000 0.00000165 0.00000018 0.00000641
11 -228.433568758067 -0.0000000000 0.00000042 0.00000005 0.00000142
*** SCF converged in 11 iterations. Time: 1.52 sec.
Spin-Restricted Hartree-Fock:
-----------------------------
Total Energy : -228.4335687581 a.u.
Electronic Energy : -387.3224182575 a.u.
Nuclear Repulsion Energy : 158.8888494994 a.u.
------------------------------------
Gradient Norm : 0.0000004225 a.u.
Ground State Information
------------------------
Charge of Molecule : 0.0
Multiplicity (2S+1) : 1.0
Magnetic Quantum Number (M_S) : 0.0
Spin Restricted Orbitals
------------------------
Molecular Orbital No. 14:
--------------------------
Occupation: 2.000 Energy: -0.58794 a.u.
( 1 C 1p+1: 0.25) ( 3 C 1p+1: -0.27) ( 4 C 1p+1: 0.30)
( 5 C 1p+1: -0.18) ( 5 C 1p-1: 0.18) ( 7 H 1s : -0.25)
( 8 H 1s : 0.25)
Molecular Orbital No. 15:
--------------------------
Occupation: 2.000 Energy: -0.56268 a.u.
( 3 C 1p-1: -0.31) ( 4 C 1p+1: 0.20) ( 4 C 1p-1: 0.24)
( 6 O 1p+1: 0.37) ( 6 O 2p+1: 0.27) ( 7 H 1s : 0.20)
( 8 H 1s : 0.20)
Molecular Orbital No. 16:
--------------------------
Occupation: 2.000 Energy: -0.54763 a.u.
( 1 C 1p+1: 0.25) ( 3 C 1p+1: -0.23) ( 3 C 1p-1: -0.20)
( 4 C 1p-1: 0.29) ( 5 C 1p+1: 0.24) ( 6 O 1p+1: -0.32)
( 6 O 2p+1: -0.24)
Molecular Orbital No. 17:
--------------------------
Occupation: 2.000 Energy: -0.39987 a.u.
( 3 C 1p0 : 0.32) ( 3 C 2p0 : 0.24) ( 4 C 1p0 : 0.32)
( 4 C 2p0 : 0.24) ( 6 O 1p0 : -0.34) ( 6 O 2p0 : -0.30)
Molecular Orbital No. 18:
--------------------------
Occupation: 2.000 Energy: -0.32333 a.u.
( 1 C 1p0 : 0.35) ( 1 C 2p0 : 0.29) ( 3 C 1p0 : 0.23)
( 3 C 2p0 : 0.21) ( 4 C 1p0 : -0.23) ( 4 C 2p0 : -0.21)
( 5 C 1p0 : -0.35) ( 5 C 2p0 : -0.29)
Molecular Orbital No. 19:
--------------------------
Occupation: 0.000 Energy: 0.14744 a.u.
( 1 C 1p0 : 0.32) ( 1 C 2p0 : 0.56) ( 3 C 1p0 : -0.18)
( 3 C 2p0 : -0.35) ( 4 C 1p0 : -0.18) ( 4 C 2p0 : -0.35)
( 5 C 1p0 : 0.32) ( 5 C 2p0 : 0.56) ( 6 O 1p0 : -0.27)
( 6 O 2p0 : -0.39)
Molecular Orbital No. 20:
--------------------------
Occupation: 0.000 Energy: 0.20445 a.u.
( 1 C 3s : 0.65) ( 1 C 2p+1: -0.30) ( 1 C 2p-1: -0.23)
( 2 H 2s : -0.98) ( 3 C 3s : 0.97) ( 3 C 2p+1: 0.29)
( 3 C 2p-1: -0.39) ( 4 C 3s : 0.97) ( 4 C 2p+1: 0.46)
( 4 C 2p-1: 0.17) ( 5 C 3s : 0.65) ( 5 C 2p-1: 0.36)
( 7 H 2s : -1.31) ( 8 H 2s : -1.31) ( 9 H 2s : -0.98)
Molecular Orbital No. 21:
--------------------------
Occupation: 0.000 Energy: 0.22624 a.u.
( 1 C 1p0 : 0.22) ( 1 C 2p0 : 0.54) ( 3 C 1p0 : -0.32)
( 3 C 2p0 : -0.83) ( 4 C 1p0 : 0.32) ( 4 C 2p0 : 0.83)
( 5 C 1p0 : -0.22) ( 5 C 2p0 : -0.54)
Molecular Orbital No. 22:
--------------------------
Occupation: 0.000 Energy: 0.23036 a.u.
( 1 C 3s : -0.73) ( 1 C 2p+1: 0.85) ( 1 C 2p-1: 0.79)
( 2 H 2s : 1.89) ( 3 C 3s : -1.11) ( 3 C 2p-1: -0.24)
( 4 C 3s : 1.11) ( 4 C 2p-1: -0.20) ( 5 C 3s : 0.73)
( 5 C 1p-1: 0.16) ( 5 C 2p+1: -0.28) ( 5 C 2p-1: 1.13)
( 7 H 2s : 0.36) ( 8 H 2s : -0.36) ( 9 H 2s : -1.89)
Molecular Orbital No. 23:
--------------------------
Occupation: 0.000 Energy: 0.23249 a.u.
( 1 C 3s : -0.93) ( 1 C 2p+1: 0.62) ( 1 C 2p-1: 0.57)
( 2 H 2s : 1.56) ( 3 C 3s : 0.60) ( 3 C 2p+1: 0.33)
( 3 C 2p-1: -0.60) ( 4 C 3s : 0.60) ( 4 C 2p+1: 0.60)
( 4 C 2p-1: 0.32) ( 5 C 3s : -0.93) ( 5 C 2p+1: 0.20)
( 5 C 2p-1: -0.81) ( 7 H 2s : -1.14) ( 8 H 2s : -1.14)
( 9 H 2s : 1.56)
Then we choose all \(\pi\) orbitals (so 5 orbitals total with 6 electrons):
space=mtp.OrbSpace(molecule,scfdrv.mol_orbs)
space.CAS(6,5)
mcscfdrv=mtp.McscfDriver()
mcscfdrv.compute(molecule,basis,space)
Active space definition:
------------------------
Number of inactive (occupied) orbitals: 15
Number of active orbitals: 5
Number of virtual orbitals: 58
This is a CASSCF wavefunction: CAS(6,5)
CI expansion:
-------------
Number of determinants: 55
MCSCF Iterations
----------------
Iter. | Average Energy | E. Change | Grad. Norm | CI Iter. | Time
---------------------------------------------------------------------
1 -228.446979427 0.0e+00 2.7e-02 1 0:00:00
2 -228.450765627 -3.8e-03 2.4e-02 1 0:00:00
3 -228.453531911 -2.8e-03 2.3e-02 1 0:00:00
4 -228.456176162 -2.6e-03 3.4e-02 1 0:00:00
5 -228.457537064 -1.4e-03 4.7e-02 1 0:00:00
6 -228.460418228 -2.9e-03 4.3e-02 1 0:00:00
7 -228.471082052 -1.1e-02 3.8e-02 1 0:00:00
8 -228.473011579 -1.9e-03 3.6e-02 1 0:00:00
9 -228.475556670 -2.5e-03 2.2e-02 1 0:00:00
10 -228.478676266 -3.1e-03 6.8e-02 1 0:00:00
11 -228.480340144 -1.7e-03 2.9e-02 1 0:00:00
12 -228.481278202 -9.4e-04 1.7e-02 1 0:00:00
13 -228.482360482 -1.1e-03 2.9e-02 1 0:00:00
14 -228.482793880 -4.3e-04 1.6e-02 1 0:00:00
15 -228.483060513 -2.7e-04 1.3e-02 1 0:00:00
16 -228.483180065 -1.2e-04 7.9e-03 1 0:00:00
17 -228.483283941 -1.0e-04 7.3e-03 1 0:00:00
18 -228.483789360 -5.1e-04 4.8e-02 1 0:00:00
19 -228.483035252 7.5e-04 1.3e-01 1 0:00:00
20 -228.484530293 -1.5e-03 9.1e-02 1 0:00:00
21 -228.486657626 -2.1e-03 6.5e-02 1 0:00:00
22 -228.489300742 -2.6e-03 5.0e-02 1 0:00:00
23 -228.490488056 -1.2e-03 2.4e-02 1 0:00:00
24 -228.490807405 -3.2e-04 1.2e-02 1 0:00:00
25 -228.490896271 -8.9e-05 7.2e-03 1 0:00:00
26 -228.490946144 -5.0e-05 4.6e-03 1 0:00:00
27 -228.490980817 -3.5e-05 4.9e-03 1 0:00:00
28 -228.491004478 -2.4e-05 4.2e-03 1 0:00:00
29 -228.491012756 -8.3e-06 2.2e-03 1 0:00:00
30 -228.491015191 -2.4e-06 1.1e-03 1 0:00:00
31 -228.491015957 -7.7e-07 7.2e-04 1 0:00:00
32 -228.491016262 -3.0e-07 4.5e-04 1 0:00:00
33 -228.491016379 -1.2e-07 2.5e-04 1 0:00:00
34 -228.491016438 -5.9e-08 2.1e-04 1 0:00:00
35 -228.491016470 -3.2e-08 1.8e-04 1 0:00:00
36 -228.491016485 -1.5e-08 1.1e-04 1 0:00:00
37 -228.491016490 -5.0e-09 5.2e-05 1 0:00:00
** Convergence reached in 37 iterations
WARNING! Your active orbitals may have changed significantly
Min_overlap = 0.05535591470376259
Final results
-------------
* State 1
- Energy: -228.49101648965615
- S^2 : 0.00 (multiplicity = 1.0 )
- Natural orbitals
1.99000 1.97270 1.97224 0.03735 0.02771
Spin Restricted Orbitals
------------------------
Molecular Orbital No. 14:
--------------------------
Occupation: 2.000 Energy: -0.54990 a.u.
( 1 C 1p+1: 0.22) ( 2 H 1s : -0.18) ( 3 C 1p+1: -0.17)
( 3 C 1p-1: -0.33) ( 3 C 2p-1: -0.15) ( 4 C 1p-1: 0.38)
( 4 C 2p-1: 0.17) ( 5 C 1p+1: 0.19)
Molecular Orbital No. 15:
--------------------------
Occupation: 2.000 Energy: -0.46675 a.u.
( 3 C 1p0 : -0.17) ( 4 C 1p0 : -0.37) ( 4 C 2p0 : -0.27)
( 5 C 1p0 : -0.36) ( 5 C 2p0 : -0.25)
Molecular Orbital No. 16:
--------------------------
Occupation: 1.990 Energy: -0.35333 a.u.
( 1 C 1p0 : 0.40) ( 1 C 2p0 : 0.30) ( 3 C 1p0 : 0.36)
( 3 C 2p0 : 0.30) ( 5 C 1p0 : -0.23) ( 5 C 2p0 : -0.22)
Molecular Orbital No. 17:
--------------------------
Occupation: 1.973 Energy: -0.81321 a.u.
( 1 C 2s : -0.17) ( 1 C 1p-1: -0.33) ( 6 O 1p+1: -0.21)
( 6 O 1p-1: 0.47) ( 6 O 2p+1: -0.15) ( 6 O 2p-1: 0.27)
Molecular Orbital No. 18:
--------------------------
Occupation: 1.972 Energy: -0.53651 a.u.
( 6 O 1p0 : -0.59) ( 6 O 2p0 : -0.43)
Molecular Orbital No. 19:
--------------------------
Occupation: 0.037 Energy: -0.02491 a.u.
( 1 C 1p0 : 0.47) ( 1 C 2p0 : 0.25) ( 3 C 1p0 : -0.29)
( 3 C 2p0 : -0.28) ( 5 C 1p0 : 0.24) ( 5 C 2p0 : 0.22)
( 6 O 1p0 : -0.66) ( 6 O 2p0 : 0.18)
Molecular Orbital No. 20:
--------------------------
Occupation: 0.028 Energy: -0.02484 a.u.
( 1 C 1s : -0.16) ( 1 C 2s : 0.49) ( 1 C 1p+1: -0.24)
( 1 C 1p-1: 0.62) ( 1 C 2p-1: 0.17) ( 6 O 2s : -0.45)
( 6 O 1p+1: -0.28) ( 6 O 1p-1: 0.73)
Molecular Orbital No. 21:
--------------------------
Occupation: 0.000 Energy: 0.20635 a.u.
( 1 C 3s : -0.74) ( 1 C 2p+1: 0.36) ( 1 C 2p-1: 0.22)
( 2 H 2s : 1.07) ( 3 C 3s : -0.95) ( 3 C 2p+1: -0.25)
( 3 C 2p-1: 0.32) ( 4 C 3s : -0.94) ( 4 C 2p+1: -0.46)
( 4 C 2p-1: -0.19) ( 5 C 3s : -0.72) ( 5 C 2p-1: -0.41)
( 7 H 2s : 1.18) ( 8 H 2s : 1.30) ( 9 H 2s : 1.08)
Molecular Orbital No. 22:
--------------------------
Occupation: 0.000 Energy: 0.22212 a.u.
( 1 C 2p0 : -0.39) ( 3 C 1p0 : 0.28) ( 3 C 2p0 : 0.72)
( 4 C 1p0 : -0.36) ( 4 C 2p0 : -0.89) ( 5 C 1p0 : 0.28)
( 5 C 2p0 : 0.65)
Molecular Orbital No. 23:
--------------------------
Occupation: 0.000 Energy: 0.23115 a.u.
( 1 C 3s : -0.57) ( 1 C 2p+1: 0.75) ( 1 C 2p-1: 0.71)
( 2 H 2s : 1.62) ( 3 C 3s : -1.21) ( 4 C 3s : 0.96)
( 4 C 2p+1: -0.27) ( 4 C 2p-1: -0.27) ( 5 C 2s : 0.15)
( 5 C 3s : 0.86) ( 5 C 1p-1: 0.17) ( 5 C 2p+1: -0.31)
( 5 C 2p-1: 1.26) ( 7 H 2s : 0.53) ( 9 H 2s : -2.12)
Total MCSCF time: 00:00:13
You see the convergence is very poor. Additionally, we get a final warning that our active space has changed significantly. The reason for it is that we started with the wrong active orbitals and the MCSCF had to make up for this. Let’s visualize the resulting active orbitals using the multipsi orbital viewer:
orbviewer=mtp.OrbitalViewer()
orbviewer.plot(molecule,basis,space)
We can see that we got 3 \(\pi\) orbitals correctly, but 2 were not. By default, when only specifying the number of electrons and orbitals like we did here, the program chooses those around the HOMO-LUMO gap. While it is sometimes fine, most often the orbitals we want are not those, and we thus need to resort to visual inspection.
Fortunately, the multipsi orbital viewer provides us with a very easy way to select the active orbitals and save them into a h5 file. Using this tool, we can remove the 2 \(\sigma\) orbitals and select the missing \(\pi\). We save the result in file “furan-cas.h5”.
space=mtp.OrbSpace(molecule,"furan-cas.h5")
mcscfdrv=mtp.McscfDriver()
mcscfdrv.compute(molecule,basis,space)
Active space definition:
------------------------
Number of inactive (occupied) orbitals: 15
Number of active orbitals: 5
Number of virtual orbitals: 58
This is a CASSCF wavefunction: CAS(6,5)
CI expansion:
-------------
Number of determinants: 55
MCSCF Iterations
----------------
Iter. | Average Energy | E. Change | Grad. Norm | CI Iter. | Time
---------------------------------------------------------------------
1 -228.469597501 0.0e+00 7.5e-02 1 0:00:00
2 -228.480395202 -1.1e-02 4.0e-02 1 0:00:00
3 -228.483481276 -3.1e-03 6.4e-02 1 0:00:00
4 -228.486710082 -3.2e-03 2.8e-02 1 0:00:00
5 -228.487821657 -1.1e-03 1.2e-02 1 0:00:00
6 -228.488001560 -1.8e-04 3.5e-03 1 0:00:00
7 -228.488013666 -1.2e-05 9.4e-04 1 0:00:00
8 -228.488014031 -3.6e-07 3.0e-04 1 0:00:00
9 -228.488014075 -4.4e-08 8.0e-05 1 0:00:00
10 -228.488014076 -1.8e-09 2.5e-05 1 0:00:00
** Convergence reached in 10 iterations
Final results
-------------
* State 1
- Energy: -228.4880140764987
- S^2 : 0.00 (multiplicity = 1.0 )
- Natural orbitals
1.99451 1.94212 1.90287 0.09707 0.06342
Spin Restricted Orbitals
------------------------
Molecular Orbital No. 14:
--------------------------
Occupation: 2.000 Energy: -0.55816 a.u.
( 3 C 1p-1: 0.32) ( 4 C 1p+1: -0.19) ( 4 C 1p-1: -0.26)
( 6 O 1p+1: -0.35) ( 6 O 2p+1: -0.26) ( 7 H 1s : -0.20)
( 8 H 1s : -0.20)
Molecular Orbital No. 15:
--------------------------
Occupation: 2.000 Energy: -0.54203 a.u.
( 1 C 1p+1: 0.25) ( 3 C 1p+1: -0.23) ( 3 C 1p-1: -0.18)
( 4 C 1p-1: 0.27) ( 5 C 1p+1: 0.24) ( 6 O 1p+1: -0.34)
( 6 O 2p+1: -0.27)
Molecular Orbital No. 16:
--------------------------
Occupation: 1.995 Energy: -0.58408 a.u.
( 6 O 1p0 : -0.59) ( 6 O 2p0 : -0.43)
Molecular Orbital No. 17:
--------------------------
Occupation: 1.942 Energy: -0.43681 a.u.
( 1 C 1p0 : 0.19) ( 3 C 1p0 : 0.34) ( 3 C 2p0 : 0.23)
( 4 C 1p0 : 0.34) ( 4 C 2p0 : 0.23) ( 5 C 1p0 : 0.19)
( 6 O 2p0 : -0.15)
Molecular Orbital No. 18:
--------------------------
Occupation: 1.903 Energy: -0.33694 a.u.
( 1 C 1p0 : 0.36) ( 1 C 2p0 : 0.28) ( 3 C 1p0 : 0.24)
( 3 C 2p0 : 0.20) ( 4 C 1p0 : -0.24) ( 4 C 2p0 : -0.20)
( 5 C 1p0 : -0.36) ( 5 C 2p0 : -0.28)
Molecular Orbital No. 19:
--------------------------
Occupation: 0.097 Energy: -0.03367 a.u.
( 1 C 1p0 : -0.48) ( 1 C 2p0 : -0.30) ( 3 C 1p0 : 0.30)
( 3 C 2p0 : 0.18) ( 4 C 1p0 : 0.30) ( 4 C 2p0 : 0.18)
( 5 C 1p0 : -0.48) ( 5 C 2p0 : -0.30) ( 6 O 1p0 : 0.28)
( 6 O 2p0 : 0.27)
Molecular Orbital No. 20:
--------------------------
Occupation: 0.063 Energy: -0.02346 a.u.
( 1 C 1p0 : 0.36) ( 1 C 2p0 : 0.24) ( 3 C 1p0 : -0.51)
( 3 C 2p0 : -0.39) ( 4 C 1p0 : 0.51) ( 4 C 2p0 : 0.39)
( 5 C 1p0 : -0.36) ( 5 C 2p0 : -0.24)
Molecular Orbital No. 21:
--------------------------
Occupation: 0.000 Energy: 0.20569 a.u.
( 1 C 3s : -0.62) ( 1 C 2p+1: 0.28) ( 1 C 2p-1: 0.21)
( 2 H 2s : 0.93) ( 3 C 3s : -0.99) ( 3 C 2p+1: -0.30)
( 3 C 2p-1: 0.41) ( 4 C 3s : -0.99) ( 4 C 2p+1: -0.48)
( 4 C 2p-1: -0.18) ( 5 C 3s : -0.62) ( 5 C 2p-1: -0.33)
( 7 H 2s : 1.35) ( 8 H 2s : 1.35) ( 9 H 2s : 0.93)
Molecular Orbital No. 22:
--------------------------
Occupation: 0.000 Energy: 0.23274 a.u.
( 1 C 3s : 0.75) ( 1 C 2p+1: -0.85) ( 1 C 2p-1: -0.79)
( 2 H 2s : -1.89) ( 3 C 3s : 1.12) ( 3 C 2p-1: 0.22)
( 4 C 3s : -1.13) ( 4 C 2p-1: 0.19) ( 5 C 3s : -0.74)
( 5 C 1p-1: -0.16) ( 5 C 2p+1: 0.28) ( 5 C 2p-1: -1.12)
( 7 H 2s : -0.40) ( 8 H 2s : 0.41) ( 9 H 2s : 1.87)
Molecular Orbital No. 23:
--------------------------
Occupation: 0.000 Energy: 0.23429 a.u.
( 1 C 3s : -0.96) ( 1 C 2p+1: 0.62) ( 1 C 2p-1: 0.58)
( 2 H 2s : 1.58) ( 3 C 3s : 0.57) ( 3 C 2p+1: 0.32)
( 3 C 2p-1: -0.59) ( 4 C 3s : 0.57) ( 4 C 2p+1: 0.59)
( 4 C 2p-1: 0.32) ( 5 C 3s : -0.96) ( 5 C 2p+1: 0.21)
( 5 C 2p-1: -0.83) ( 7 H 2s : -1.10) ( 8 H 2s : -1.10)
( 9 H 2s : 1.60)
Total MCSCF time: 00:00:04
The calculation converged much faster, and you can also notice that the occupation numbers are a bit further away from 2 and 0, which is common of \(\pi\) bonds. Let us visualize them to confirm:
orbviewer=mtp.OrbitalViewer()
orbviewer.plot(molecule,basis,space)
It is always recommended to inspect the orbitals after your calculations, as the optimisation may sometimes swap our active orbitals, especially if their occupations are close to 2 or 0.
State-specific and state-averaged CASSCF¶
When one is interested in more than just the ground state, for example in spectroscopy or photochemistry, we can specify this in the mcscf driver:
nstates=3
mcscfdrv.compute(molecule,basis,space, nstates)
Active space definition:
------------------------
Number of inactive (occupied) orbitals: 15
Number of active orbitals: 5
Number of virtual orbitals: 58
This is a CASSCF wavefunction: CAS(6,5)
CI expansion:
-------------
Number of determinants: 55
MCSCF Iterations
----------------
Iter. | Average Energy | E. Change | Grad. Norm | CI Iter. | Time
---------------------------------------------------------------------
1 -228.278942411 0.0e+00 1.3e-01 1 0:00:00
2 -228.296679706 -1.8e-02 5.9e-02 1 0:00:00
3 -228.298696659 -2.0e-03 2.1e-02 1 0:00:00
4 -228.298795467 -9.9e-05 5.2e-03 1 0:00:00
5 -228.298804959 -9.5e-06 1.1e-03 1 0:00:00
6 -228.298805381 -4.2e-07 2.5e-04 1 0:00:00
7 -228.298805417 -3.6e-08 7.2e-05 1 0:00:00
8 -228.298805419 -2.7e-09 2.2e-05 1 0:00:00
** Convergence reached in 8 iterations
Final results
-------------
* State 1
- Energy: -228.48189324195783
- S^2 : 0.00 (multiplicity = 1.0 )
- Natural orbitals
1.99296 1.94507 1.90727 0.09336 0.06134
* State 2
- Energy: -228.23493686018864
- S^2 : -0.00 (multiplicity = 1.0 )
- Natural orbitals
1.99023 1.43277 1.34328 0.74246 0.49127
* State 3
- Energy: -228.179586155675
- S^2 : 0.00 (multiplicity = 1.0 )
- Natural orbitals
1.97962 1.95162 1.02709 0.97835 0.06332
Spin Restricted Orbitals
------------------------
Molecular Orbital No. 14:
--------------------------
Occupation: 2.000 Energy: -0.56289 a.u.
( 3 C 1p-1: 0.30) ( 4 C 1p+1: -0.21) ( 4 C 1p-1: -0.22)
( 6 O 1p+1: -0.39) ( 6 O 2p+1: -0.29) ( 7 H 1s : -0.20)
( 8 H 1s : -0.20)
Molecular Orbital No. 15:
--------------------------
Occupation: 2.000 Energy: -0.54801 a.u.
( 1 C 1p+1: 0.25) ( 2 H 1s : -0.15) ( 3 C 1p+1: -0.22)
( 3 C 1p-1: -0.22) ( 4 C 1p-1: 0.31) ( 5 C 1p+1: 0.23)
( 6 O 1p+1: -0.29) ( 6 O 2p+1: -0.22) ( 9 H 1s : -0.15)
Molecular Orbital No. 16:
--------------------------
Occupation: 1.975 Energy: -0.60369 a.u.
( 1 C 1p0 : -0.18) ( 5 C 1p0 : -0.18) ( 6 O 1p0 : -0.54)
( 6 O 2p0 : -0.38)
Molecular Orbital No. 17:
--------------------------
Occupation: 1.762 Energy: -0.38033 a.u.
( 3 C 1p0 : -0.36) ( 3 C 2p0 : -0.26) ( 4 C 1p0 : -0.36)
( 4 C 2p0 : -0.26) ( 6 O 1p0 : 0.23) ( 6 O 2p0 : 0.20)
Molecular Orbital No. 18:
--------------------------
Occupation: 1.412 Energy: -0.24773 a.u.
( 1 C 1p0 : 0.37) ( 1 C 2p0 : 0.30) ( 3 C 1p0 : 0.22)
( 3 C 2p0 : 0.18) ( 4 C 1p0 : -0.22) ( 4 C 2p0 : -0.18)
( 5 C 1p0 : -0.37) ( 5 C 2p0 : -0.30)
Molecular Orbital No. 19:
--------------------------
Occupation: 0.632 Energy: -0.05471 a.u.
( 1 C 1p0 : -0.38) ( 1 C 2p0 : -0.48) ( 3 C 1p0 : 0.18)
( 3 C 2p0 : 0.28) ( 4 C 1p0 : 0.18) ( 4 C 2p0 : 0.28)
( 5 C 1p0 : -0.38) ( 5 C 2p0 : -0.48) ( 6 O 1p0 : 0.33)
( 6 O 2p0 : 0.37)
Molecular Orbital No. 20:
--------------------------
Occupation: 0.220 Energy: -0.03070 a.u.
( 1 C 1p0 : -0.28) ( 1 C 2p0 : -0.38) ( 3 C 1p0 : 0.45)
( 3 C 2p0 : 0.58) ( 4 C 1p0 : -0.45) ( 4 C 2p0 : -0.58)
( 5 C 1p0 : 0.28) ( 5 C 2p0 : 0.38)
Molecular Orbital No. 21:
--------------------------
Occupation: 0.000 Energy: 0.20516 a.u.
( 1 C 3s : -0.63) ( 1 C 2p+1: 0.29) ( 1 C 2p-1: 0.23)
( 2 H 2s : 0.95) ( 3 C 3s : -0.98) ( 3 C 2p+1: -0.30)
( 3 C 2p-1: 0.40) ( 4 C 3s : -0.98) ( 4 C 2p+1: -0.47)
( 4 C 2p-1: -0.17) ( 5 C 3s : -0.64) ( 5 C 2p-1: -0.35)
( 7 H 2s : 1.33) ( 8 H 2s : 1.33) ( 9 H 2s : 0.95)
Molecular Orbital No. 22:
--------------------------
Occupation: 0.000 Energy: 0.23176 a.u.
( 1 C 3s : 0.74) ( 1 C 2p+1: -0.85) ( 1 C 2p-1: -0.79)
( 2 H 2s : -1.88) ( 3 C 3s : 1.11) ( 3 C 2p-1: 0.21)
( 4 C 3s : -1.10) ( 4 C 2p-1: 0.19) ( 5 C 3s : -0.75)
( 5 C 1p-1: -0.16) ( 5 C 2p+1: 0.28) ( 5 C 2p-1: -1.13)
( 7 H 2s : -0.39) ( 8 H 2s : 0.38) ( 9 H 2s : 1.89)
Molecular Orbital No. 23:
--------------------------
Occupation: 0.000 Energy: 0.23380 a.u.
( 1 C 3s : 0.96) ( 1 C 2p+1: -0.62) ( 1 C 2p-1: -0.57)
( 2 H 2s : -1.58) ( 3 C 3s : -0.58) ( 3 C 2p+1: -0.32)
( 3 C 2p-1: 0.59) ( 4 C 3s : -0.59) ( 4 C 2p+1: -0.59)
( 4 C 2p-1: -0.32) ( 5 C 3s : 0.96) ( 5 C 2p+1: -0.21)
( 5 C 2p-1: 0.81) ( 7 H 2s : 1.11) ( 8 H 2s : 1.11)
( 9 H 2s : -1.57)
Total MCSCF time: 00:00:03
By default, the program then uses a state-averaged optimisation. This means the orbitals are found to minimize the average of all the selected states (here 3). By doing so, we obtain a compromise and all states are treated at a similar accuracy, favouring energy cancellation in the energy differences (i.e. excitation energies). On the other hand, the orbitals are not optimal for any specific state, and here we can see for example that the ground state energy is slightly higher than the one found in the previous section.
Alternatively, one can decide to optimize the orbitals for a specific state, or biasing the average for a specific state. For this, we use the state-average optimizer but with custom weights instead of the default equal weights:
mcscfdrv.compute(molecule,basis,space, nstates, [0,1,0])
Active space definition:
------------------------
Number of inactive (occupied) orbitals: 15
Number of active orbitals: 5
Number of virtual orbitals: 58
This is a CASSCF wavefunction: CAS(6,5)
CI expansion:
-------------
Number of determinants: 55
MCSCF Iterations
----------------
Iter. | Average Energy | E. Change | Grad. Norm | CI Iter. | Time
---------------------------------------------------------------------
1 -228.234936860 0.0e+00 1.1e-01 1 0:00:00
2 -228.239446480 -4.5e-03 3.8e-02 1 0:00:00
3 -228.240111316 -6.6e-04 1.1e-02 1 0:00:00
4 -228.240166877 -5.6e-05 4.0e-03 1 0:00:00
5 -228.240174438 -7.6e-06 9.4e-04 1 0:00:00
6 -228.240174649 -2.1e-07 2.7e-04 1 0:00:00
7 -228.240174664 -1.5e-08 3.9e-05 1 0:00:00
8 -228.240174665 -5.6e-10 1.1e-05 1 0:00:00
** Convergence reached in 8 iterations
Final results
-------------
* State 1
- Energy: -228.47917835455752
- S^2 : -0.00 (multiplicity = 1.0 )
- Natural orbitals
1.99448 1.94016 1.89395 0.10619 0.06522
* State 2
- Energy: -228.24017466465176
- S^2 : -0.00 (multiplicity = 1.0 )
- Natural orbitals
1.99260 1.43909 1.31021 0.75252 0.50558
* State 3
- Energy: -228.16190451234334
- S^2 : 0.00 (multiplicity = 1.0 )
- Natural orbitals
1.97985 1.97985 1.02838 0.97599 0.06357
Spin Restricted Orbitals
------------------------
Molecular Orbital No. 14:
--------------------------
Occupation: 2.000 Energy: -0.56225 a.u.
( 3 C 1p-1: 0.35) ( 4 C 1p+1: -0.16) ( 4 C 1p-1: -0.31)
( 6 O 1p+1: -0.28) ( 6 O 2p+1: -0.20) ( 7 H 1s : -0.19)
( 8 H 1s : -0.19)
Molecular Orbital No. 15:
--------------------------
Occupation: 2.000 Energy: -0.54175 a.u.
( 1 C 1p+1: 0.23) ( 3 C 1p+1: -0.22) ( 4 C 1p-1: 0.21)
( 5 C 1p+1: 0.23) ( 6 O 1p+1: -0.41) ( 6 O 2p+1: -0.31)
Molecular Orbital No. 16:
--------------------------
Occupation: 1.993 Energy: -0.59375 a.u.
( 6 O 1p0 : -0.56) ( 6 O 2p0 : -0.42)
Molecular Orbital No. 17:
--------------------------
Occupation: 1.439 Energy: -0.32504 a.u.
( 3 C 1p0 : -0.40) ( 3 C 2p0 : -0.31) ( 4 C 1p0 : -0.40)
( 4 C 2p0 : -0.31) ( 6 O 1p0 : 0.15)
Molecular Orbital No. 18:
--------------------------
Occupation: 1.310 Energy: -0.25299 a.u.
( 1 C 1p0 : 0.43) ( 1 C 2p0 : 0.36) ( 5 C 1p0 : -0.43)
( 5 C 2p0 : -0.36)
Molecular Orbital No. 19:
--------------------------
Occupation: 0.753 Energy: -0.10340 a.u.
( 1 C 1p0 : -0.46) ( 1 C 2p0 : -0.43) ( 3 C 2p0 : 0.21)
( 4 C 2p0 : 0.21) ( 5 C 1p0 : -0.46) ( 5 C 2p0 : -0.43)
( 6 O 1p0 : 0.30) ( 6 O 2p0 : 0.37)
Molecular Orbital No. 20:
--------------------------
Occupation: 0.506 Energy: -0.05610 a.u.
( 1 C 1p0 : 0.18) ( 1 C 2p0 : 0.34) ( 3 C 1p0 : -0.48)
( 3 C 2p0 : -0.63) ( 4 C 1p0 : 0.48) ( 4 C 2p0 : 0.63)
( 5 C 1p0 : -0.18) ( 5 C 2p0 : -0.34)
Molecular Orbital No. 21:
--------------------------
Occupation: 0.000 Energy: 0.20391 a.u.
( 1 C 3s : -0.51) ( 1 C 2p+1: 0.21) ( 1 C 2p-1: 0.15)
( 2 H 2s : 0.74) ( 3 C 3s : -1.05) ( 3 C 2p+1: -0.34)
( 3 C 2p-1: 0.47) ( 4 C 3s : -1.05) ( 4 C 2p+1: -0.54)
( 4 C 2p-1: -0.21) ( 5 C 3s : -0.51) ( 5 C 2p-1: -0.24)
( 7 H 2s : 1.46) ( 8 H 2s : 1.46) ( 9 H 2s : 0.74)
Molecular Orbital No. 22:
--------------------------
Occupation: 0.000 Energy: 0.23646 a.u.
( 1 C 3s : 0.84) ( 1 C 2p+1: -0.73) ( 1 C 2p-1: -0.64)
( 2 H 2s : -1.71) ( 3 C 3s : 1.16) ( 3 C 2p+1: 0.28)
( 4 C 3s : -1.17) ( 4 C 2p+1: -0.28) ( 5 C 3s : -0.83)
( 5 C 2p+1: 0.26) ( 5 C 2p-1: -0.93) ( 7 H 2s : -0.89)
( 8 H 2s : 0.89) ( 9 H 2s : 1.69)
Molecular Orbital No. 23:
--------------------------
Occupation: 0.000 Energy: 0.23718 a.u.
( 1 C 3s : 1.05) ( 1 C 2p+1: -0.64) ( 1 C 2p-1: -0.60)
( 2 H 2s : -1.68) ( 3 C 3s : -0.47) ( 3 C 2p+1: -0.27)
( 3 C 2p-1: 0.54) ( 4 C 3s : -0.46) ( 4 C 2p+1: -0.52)
( 4 C 2p-1: -0.30) ( 5 C 3s : 1.06) ( 5 C 2p+1: -0.21)
( 5 C 2p-1: 0.86) ( 7 H 2s : 0.94) ( 8 H 2s : 0.93)
( 9 H 2s : -1.70)
Total MCSCF time: 00:00:03
Here we selected 0 weights on all states except the second one, and as expected, that state is now lower in energy than in the state-averaged case.
Computing each state this way (called the state-specific optimization) provides the best description of each state but suffers from 2 major issues. First, it is not always easy to optimize excited states this way, as the order of the states may change depending on the orbitals and it is thus difficult to truly follow a specific state. Additionally, the orthonormality of the states wavefunction is not enforced, and thus the excited states may appear lower in energy than they should be. This can be fixed a posteriori by computing the overlap and Hamiltonian between the different states and